function single_level_analysis_visualization_functions(datain, visid, opt, flag)
    %%% global parameters %%%
    Fsi = 20;
    Fsb = 20;
    tstep = [2, 10];
    crg = [0, 0.002];
    clrs = [0, 0.4470, 0.7410; 0.8500, 0.3250, 0.0980; 0.9290, 0.6940, 0.1250; 0.4940, 0.1840, 0.5560; 0.4660, 0.6740, 0.1880; 0.3010, 0.7450, 0.9330; 0.6350, 0.0780, 0.1840; 0.25, 0.25, 0.25];
%     clrs = [0, 0.4470, 0.7410; 0.8500, 0.3250, 0.0980; 0, 0.4470, 0.7410; 0.8500, 0.3250, 0.0980; 0, 0.4470, 0.7410; 0.8500, 0.3250, 0.0980;];

    switch visid
        case 'all_neuron_trial'
            %%% local parameters %%%
            datas = datain{1};
            [nc, nw, nm] = size(datas);
            
            %%% main %%%
            figure(gcf)
            clf
            mn = zeros(nc, nw, nm);
            mx = zeros(nc, nw, nm);
            h = cell(nc, nw, nm);
            for i = 1: nc
                for j = 1: nw
                    for k = 1: nm
                        h{i, j, k} = subplot(nc * nw, nm, k + (i - 1) * nw * nm + (j - 1) * nm);
                        tp = datas{i, j, k}{flag};
                        if ~isempty(tp)
                            ttmp = squeeze(mean(tp, 1))';
                            if isvector(ttmp)
                                ttmp = ttmp';
                            end
                            ttmp = ttmp - mean(ttmp(:, 1: tstep(1) * Fsi), 2);
                            [~, idx] = max(ttmp, [], 2);
                            [~, idx] = sort(idx);
                            %                         mn(i, j, k) = min(ttmp(:));
                            %                         mx(i, j, k) = max(ttmp(:));
                            mn(i, j, k) = prctile(ttmp(:), 1);
                            mx(i, j, k) = prctile(ttmp(:), 99);
                            imagesc([-tstep(1), tstep(2)], [1, length(idx)], ttmp(idx, :))
                            hold on
                            rgy = get(gca, 'ylim');
                            plot([0, 0], rgy, 'r')
                            hold off
                        end
                    end
                end
            end
            colormap('jet')
            for k = 1: nm
                mnc = min(min(mn(:, :, k)));
                mxc = max(max(mx(:, :, k)));
                for i = 1: nc
                    for j = 1: nw
%                         caxis(h{i, j, k}, [mnc, mxc])
                        caxis(h{i, j, k}, crg)
                    end
                end
            end
            
        case 'average_trial_all_neuron'
            %%% local parameters %%%
            datas = datain{1};
            [nc, nw, nm] = size(datas);
            
            %%% main %%%
            figure(gcf)
            clf
            mn = zeros(nc, nw, nm);
            mx = zeros(nc, nw, nm);
            h = cell(nc, nw, nm);
            for i = 1: nc
                for j = 1: nw
                    for k = 1: nm
                        h{i, j, k} = subplot(nc * nw, nm, k + (i - 1) * nw * nm + (j - 1) * nm);
                        tp = datas{i, j, k}{flag};
                        if ~isempty(tp)
                            ttmp = squeeze(mean(tp, 3));
                            if isvector(ttmp)
                                ttmp = ttmp';
                            end
                            ttmp = ttmp - mean(ttmp(:, 1: tstep(1) * Fsi), 2);
                            [~, idx] = max(ttmp, [], 2);
                            [~, idx] = sort(idx);
                            %                         mn(i, j, k) = min(ttmp(:));
                            %                         mx(i, j, k) = max(ttmp(:));
                            mn(i, j, k) = prctile(ttmp(:), 1);
                            mx(i, j, k) = prctile(ttmp(:), 99);
                            imagesc([-tstep(1), tstep(2)], [1, length(idx)], ttmp(idx, :))
                            hold on
                            rgy = get(gca, 'ylim');
                            plot([0, 0], rgy, 'r')
                            hold off
                        end
                    end
                end
            end
            colormap('jet')
            for k = 1: nm
                mnc = min(min(mn(:, :, k)));
                mxc = max(max(mx(:, :, k)));
                for i = 1: nc
                    for j = 1: nw
%                         caxis(h{i, j, k}, [mnc, mxc])
%                         caxis(h{i, j, k}, [0, mxc / 6])
                        caxis(h{i, j, k}, crg)
                    end
                end
            end
                        
        case 'all_trial_all_neuron'
            %%% local parameters %%%
            datas = datain{1};
            [nc, nw, nm] = size(datas);
%             crg = [0, 0.002];
            
            %%% main %%%
            figure(gcf)
            clf
            mn = zeros(nc, nw, nm);
            mx = zeros(nc, nw, nm);
            h = cell(nc, nw, nm);
            for i = 1: nc
                for j = 1: nw
                    for k = 1: nm
                        h{i, j, k} = subplot(nc * nw, nm, k + (i - 1) * nw * nm + (j - 1) * nm);
                        tp = datas{i, j, k}{flag};
                        if ~isempty(tp)
                            ttmp = reshape(permute(tp, [2, 1, 3]), size(tp, 2), [])';
                            if isvector(ttmp)
                                ttmp = ttmp';
                            end
                            ttmp = ttmp - mean(ttmp(:, 1: tstep(1) * Fsi), 2);
                            [~, idx] = max(ttmp, [], 2);
                            [~, idx] = sort(idx);
                            %                         mn(i, j, k) = min(ttmp(:));
                            %                         mx(i, j, k) = max(ttmp(:));
                            mn(i, j, k) = prctile(ttmp(:), 1);
                            mx(i, j, k) = prctile(ttmp(:), 99);
                            imagesc([-tstep(1), tstep(2)], [1, length(idx)], ttmp(idx, :))
                            hold on
                            rgy = get(gca, 'ylim');
                            plot([0, 0], rgy, 'r')
                            hold off
                        end
                    end
                end
            end
            colormap('jet')
            for k = 1: nm
                mnc = min(min(mn(:, :, k)));
                mxc = max(max(mx(:, :, k)));
                for i = 1: nc
                    for j = 1: nw
%                         caxis(h{i, j, k}, [mnc, mxc])
                        caxis(h{i, j, k}, crg)
                    end
                end
            end
            
        case 'whisk_no_whisk_trial'
            %%% local parameters %%%
%             crg = [0, 0.002];
            dataswhnwh = datain{1};
            [nc, nw, nm] = size(dataswhnwh);
            nwc = length(dataswhnwh{1, 1, 1}{flag});
            
            %%% main %%%
            figure(gcf)
            clf
            mn = zeros(nc, nw, nm, nwc);
            mx = zeros(nc, nw, nm, nwc);
            h = cell(nc, nw, nm, nwc);
            for i = 1: nc
                for j = 1: nw
                    for k = 1: nm
                        for kk = 1: nwc
                            h{i, j, k, kk} = subplot(nc * nw, nm * 3, (k - 1) * 3 + kk + ((i - 1) * nw + j - 1) * nm * 3);
                            tp = dataswhnwh{i, j, k}{flag}{kk};
                            ttmp = squeeze(mean(tp, 1))';
                            if isvector(ttmp)
                                ttmp = ttmp(:)';
                            end
                            ttmp = ttmp - mean(ttmp(:, 1: tstep(1) * Fsi), 2);
                            [~, idx] = max(ttmp, [], 2);
                            [~, idx] = sort(idx);
                            mn(i, j, k, kk) = prctile(ttmp(:), 1);
                            mx(i, j, k, kk) = prctile(ttmp(:), 99);
                            imagesc([-tstep(1), tstep(2)], [1, length(idx)], ttmp(idx, :))
                            hold on
                            rgy = get(gca, 'ylim');
                            plot([0, 0], rgy, 'r')
                            hold off
                        end
                    end
                end
            end
            colormap('jet')
            for k = 1: nm
                mnc = min(min(min(squeeze(mn(:, :, k, :)))));
                mxc = max(max(max(squeeze(mx(:, :, k, :)))));
                for i = 1: nc
                    for j = 1: nw
                        for kk = 1: nwc
%                             caxis(h{i, j, k, kk}, [mnc, mxc])
                            caxis(h{i, j, k, kk}, crg)
                        end
                    end
                end
            end
                        
        case 'whisk_no_whisk_neuron'
%             crg = [0, 0.002];
            %%% local parameters %%%
            dataswhnwh = datain{1};
            [nc, nw, nm] = size(dataswhnwh);
            nwc = length(dataswhnwh{1, 1, 1}{flag});
            
            %%% main %%%
            figure(gcf)
            clf
            mn = zeros(nc, nw, nm, nwc);
            mx = zeros(nc, nw, nm, nwc);
            h = cell(nc, nw, nm, nwc);
            for i = 1: nc
                for j = 1: nw
                    for k = 1: nm
                        for kk = 1: nwc
                            h{i, j, k, kk} = subplot(nc * nw, nm * 3, (k - 1) * 3 + kk + ((i - 1) * nw + j - 1) * nm * 3);
                            tp = dataswhnwh{i, j, k}{flag}{kk};
                            ttmp = squeeze(mean(tp, 3));
                            if isvector(ttmp)
                                ttmp = ttmp(:)';
                            end
                            ttmp = ttmp - mean(ttmp(:, 1: tstep(1) * Fsi), 2);
                            idx = 1: size(ttmp, 1);
%                             [~, idx] = max(ttmp, [], 2);
%                             [~, idx] = sort(idx);
                            mn(i, j, k, kk) = prctile(ttmp(:), 1);
                            mx(i, j, k, kk) = prctile(ttmp(:), 99);
                            imagesc([-tstep(1), tstep(2)], [1, length(idx)], ttmp(idx, :))
                            hold on
                            rgy = get(gca, 'ylim');
                            plot([0, 0], rgy, 'r')
                            hold off
                        end
                    end
                end
            end
            colormap('jet')
            for k = 1: nm
                mnc = min(min(min(squeeze(mn(:, :, k, :)))));
                mxc = max(max(max(squeeze(mx(:, :, k, :)))));
                for i = 1: nc
                    for j = 1: nw
                        for kk = 1: nwc
%                             caxis(h{i, j, k, kk}, [mnc, mxc])
                            caxis(h{i, j, k, kk}, crg)
                        end
                    end
                end
            end
            
        case 'whisk_no_whisk_average'
            %%% local parameters %%%
            dataswhnwh = datain{1};
            [nc, nw, nm] = size(dataswhnwh);
            nwc = length(dataswhnwh{1, 1, 1}{flag});
            
            %%% main %%%
            figure(gcf)
            clf
            mn = zeros(nc, nw, nm, nwc);
            mx = zeros(nc, nw, nm, nwc);
            h = cell(nc, nw, nm);
            for i = 1: nc
                for j = 1: nw
                    for k = 1: nm
                        h{i, j, k} = subplot(nc * nw, nm, k + ((i - 1) * nw + j - 1) * nm);
                        cls = {'k', 'r'};
                        hold on
                        for kk = 1: nwc
                            tp = dataswhnwh{i, j, k}{flag}{kk};
                            ttmp = squeeze(mean(mean(tp, 3), 1));
                            if isvector(ttmp)
                                ttmp = ttmp(:)';
                            end
                            ttmp = ttmp - mean(ttmp(:, 1: tstep(1) * Fsi), 2);
                            mn(i, j, k, kk) = prctile(ttmp(:), 1);
                            mx(i, j, k, kk) = prctile(ttmp(:), 99);
                            plot(-tstep(1): (1 / Fsi): tstep(2), ttmp, cls{kk}, 'linewidth', 2)
                        end
                    end
                end
            end
            colormap('jet')
            for k = 1: nm
                mnc = min(min(min(squeeze(mn(:, :, k, :)))));
                mxc = max(max(max(squeeze(mx(:, :, k, :)))));
                for i = 1: nc
                    for j = 1: nw
                        for kk = 1: nwc
                            xlim(h{i, j, k}, [-tstep(1), tstep(2)])
                            ylim(h{i, j, k}, [mnc, mxc])
                            rgy = get(h{i, j, k}, 'ylim');
                            plot(h{i, j, k}, [0, 0], rgy, 'r')
%                             hold off
                        end
                    end
                end
            end
            
        case 'wipe_no_wipe_trial'
            %%% local parameters %%%
            dataswpnwp = datain{1};
            [nc, nw, nm] = size(dataswpnwp);
            nwc = length(dataswpnwp{1, 1, 1}{flag});
            
            %%% main %%%
            figure(gcf)
            clf
            mn = zeros(nc, nw, nm, nwc);
            mx = zeros(nc, nw, nm, nwc);
            h = cell(nc, nw, nm, nwc);
            for i = 1: nc
                for j = 1: nw
                    for k = 1: nm
                        for kk = 1: nwc
                            h{i, j, k, kk} = subplot(nc * nw, nm * 3, (k - 1) * 3 + kk + ((i - 1) * nw + j - 1) * nm * 3);
                            tp = dataswpnwp{i, j, k}{flag}{kk};
                            ttmp = squeeze(mean(tp, 1))';
                            if isvector(ttmp)
                                ttmp = ttmp(:)';
                            end
                            ttmp = ttmp - mean(ttmp(:, 1: tstep(1) * Fsi), 2);
                            [~, idx] = max(ttmp, [], 2);
                            [~, idx] = sort(idx);
                            mn(i, j, k, kk) = prctile(ttmp(:), 1);
                            mx(i, j, k, kk) = prctile(ttmp(:), 99);
                            imagesc([-tstep(1), tstep(2)], [1, length(idx)], ttmp(idx, :))
                            hold on
                            rgy = get(gca, 'ylim');
                            plot([0, 0], rgy, 'r')
                            hold off
                        end
                    end
                end
            end
            colormap('jet')
            for k = 1: nm
                mnc = min(min(min(squeeze(mn(:, :, k, :)))));
                mxc = max(max(max(squeeze(mx(:, :, k, :)))));
                for i = 1: nc
                    for j = 1: nw
                        for kk = 1: nwc
                            caxis(h{i, j, k, kk}, [mnc, mxc])
                        end
                    end
                end
            end
            
        case 'wipe_no_wipe_neuron'
            %%% local parameters %%%
            dataswpnwp = datain{1};
            [nc, nw, nm] = size(dataswpnwp);
            nwc = length(dataswpnwp{1, 1, 1}{flag});
            
            %%% main %%%
            figure(gcf)
            clf
            mn = zeros(nc, nw, nm, nwc);
            mx = zeros(nc, nw, nm, nwc);
            h = cell(nc, nw, nm, nwc);
            for i = 1: nc
                for j = 1: nw
                    for k = 1: nm
                        for kk = 1: nwc
                            h{i, j, k, kk} = subplot(nc * nw, nm * 3, (k - 1) * 3 + kk + ((i - 1) * nw + j - 1) * nm * 3);
                            tp = dataswpnwp{i, j, k}{flag}{kk};
                            ttmp = squeeze(mean(tp, 3));
                            if isvector(ttmp)
                                ttmp = ttmp(:)';
                            end
                            ttmp = ttmp - mean(ttmp(:, 1: tstep(1) * Fsi), 2);
                            [~, idx] = max(ttmp, [], 2);
                            [~, idx] = sort(idx);
                            mn(i, j, k, kk) = prctile(ttmp(:), 1);
                            mx(i, j, k, kk) = prctile(ttmp(:), 99);
                            imagesc([-tstep(1), tstep(2)], [1, length(idx)], ttmp(idx, :))
                            hold on
                            rgy = get(gca, 'ylim');
                            plot([0, 0], rgy, 'r')
                            hold off
                        end
                    end
                end
            end
            colormap('jet')
            for k = 1: nm
                mnc = min(min(min(squeeze(mn(:, :, k, :)))));
                mxc = max(max(max(squeeze(mx(:, :, k, :)))));
                for i = 1: nc
                    for j = 1: nw
                        for kk = 1: nwc
                            caxis(h{i, j, k, kk}, [mnc, mxc])
                        end
                    end
                end
            end
            
        case 'wipe_no_wipe_average'
            %%% local parameters %%%
            dataswpnwp = datain{1};
            [nc, nw, nm] = size(dataswpnwp);
            nwc = length(dataswpnwp{1, 1, 1}{flag});
            
            %%% main %%%
            figure(gcf)
            clf
            mn = zeros(nc, nw, nm, nwc);
            mx = zeros(nc, nw, nm, nwc);
            h = cell(nc, nw, nm);
            for i = 1: nc
                for j = 1: nw
                    for k = 1: nm
                        h{i, j, k} = subplot(nc * nw, nm, k + ((i - 1) * nw + j - 1) * nm);
                        cls = {'k', 'r'};
                        hold on
                        for kk = 1: nwc
                            tp = dataswpnwp{i, j, k}{flag}{kk};
                            ttmp = squeeze(mean(mean(tp, 3), 1));
                            if isvector(ttmp)
                                ttmp = ttmp(:)';
                            end
                            ttmp = ttmp - mean(ttmp(:, 1: tstep(1) * Fsi), 2);
                            mn(i, j, k, kk) = prctile(ttmp(:), 1);
                            mx(i, j, k, kk) = prctile(ttmp(:), 99);
                            plot(-tstep(1): (1 / Fsi): tstep(2), ttmp, cls{kk}, 'linewidth', 2)
                        end
                    end
                end
            end
            colormap('jet')
            for k = 1: nm
                mnc = min(min(min(squeeze(mn(:, :, k, :)))));
                mxc = max(max(max(squeeze(mx(:, :, k, :)))));
                for i = 1: nc
                    for j = 1: nw
                        for kk = 1: nwc
                            xlim(h{i, j, k}, [-tstep(1), tstep(2)])
                            ylim(h{i, j, k}, [mnc, mxc])
                            rgy = get(h{i, j, k}, 'ylim');
                            plot(h{i, j, k}, [0, 0], rgy, 'r')
%                             hold off
                        end
                    end
                end
            end
            
        case 'whisk_no_whisk_mice_combined_with_stats'
            %%% local parameters %%%
            dataswhnwhmc = datain{1};
            dataswhnwhwpsmc = datain{2};
            [nc, nw, nwc] = size(dataswhnwhmc);
            scl = [1, 0.1, 15];
%             ylims = [-0.15, 0.15] * scl(flag);
%             ylims = [0, 0.15] * scl(flag);
            ylims = [0, 1] * scl(flag);
            
            %%% main %%%
            figure(gcf)
            clf
            for i = 1: nc
                for j = 1: nw
                    for k = 1: nwc
                        dtc = dataswhnwhmc{i, j, k}{flag};
                        for ii = 1: size(dtc, 1)
                            dtc(ii, :) = normalize(dtc(ii, :));
                        end
                        dtc = dtc - mean(dtc(:, 1: 2 * Fsi), 2);
                        wptc = dataswhnwhwpsmc{i, j, k};
%                         [~, id] = max(dtc, [], 2);
                        id = max(dtc, [], 2);
%                         [~, id] = max(dtc > 0, [], 2);
                        [~, id] = sort(id);
                        dtc = dtc(id, :);

%                         subplot(4 * nc, nw * nwc, (j - 1) * nwc + k + 4 * (i - 1) * nw * nwc)
%                         hold on
%                         for kk = 1: length(wptc)
%                             plot([wptc(kk), wptc(kk)] / Fsi - tstep(1), [kk - 1, kk], 'color', [0.4, 0.6, 0.8], 'linewidth', 2)
%                         end
%                         hold off
%                         xlim([-tstep(1), tstep(2)])
%                         ylim([0, length(wptc)])

                        subplot(4 * nc, nw * nwc, (j - 1) * nwc + k + 4 * (i - 1) * nw * nwc + nw * nwc * [1, 2, 3])
%                         if i == 2
%                             imagesc([-1, 2], [1, size(dtc, 1)], dtc(:, (tstep(1) - 1) * Fsi + 1: (2 + tstep(1)) * Fsi) .^ 1, ylims)
%                         else
                            imagesc([-tstep(1), tstep(2)], [1, size(dtc, 1)], dtc .^ 1, ylims)
%                             imagesc([-1, 5], [1, size(dtc, 1)], dtc(:, (tstep(1) - 1) * Fsi + 1: (5 + tstep(1)) * Fsi) .^ 1, ylims)
%                         end
                        hold on
                        plot([0, 0], get(gca, 'ylim'), 'r', 'linewidth', 1)
                        hold off
                        colormap('jet')
                    end
                end
            end
            
        case 'whisk_no_whisk_mice_combined_with_stats_mean'
            %%% local parameters %%%
            dataswhnwhmc = datain{1};
            dataswhnwhwpsmc = datain{2};
            [nc, nw, nwc] = size(dataswhnwhmc);
            scl = [1, 0.1, 15];
%             ylims = [-0.15, 0.15] * scl(flag);
%             ylims = [0, 0.15] * scl(flag);
            ylims = [0, 1] * scl(flag);
            
            %%% main %%%
            figure(gcf)
            clf
            for i = 1: nc
                for j = 1: nw
                    for k = 1: nwc
                        dtc = dataswhnwhmc{i, j, k}{flag};
                        for ii = 1: size(dtc, 1)
                            dtc(ii, :) = normalize(dtc(ii, :));
                        end
                        dtc = dtc - mean(dtc(:, 1: 2 * Fsi), 2);
                        wptc = dataswhnwhwpsmc{i, j, k};
%                         [~, id] = max(dtc, [], 2);
                        id = max(dtc, [], 2);
%                         [~, id] = max(dtc > 0, [], 2);
                        [~, id] = sort(id);
                        dtc = dtc(id, :);

%                         subplot(4 * nc, nw * nwc, (j - 1) * nwc + k + 4 * (i - 1) * nw * nwc)
%                         hold on
%                         for kk = 1: length(wptc)
%                             plot([wptc(kk), wptc(kk)] / Fsi - tstep(1), [kk - 1, kk], 'color', [0.4, 0.6, 0.8], 'linewidth', 2)
%                         end
%                         hold off
%                         xlim([-tstep(1), tstep(2)])
%                         ylim([0, length(wptc)])

                        subplot(4 * nc, nw * nwc, (j - 1) * nwc + k + 4 * (i - 1) * nw * nwc + nw * nwc * [1, 2, 3])
%                         if i == 2
%                             imagesc([-1, 2], [1, size(dtc, 1)], dtc(:, (tstep(1) - 1) * Fsi + 1: (2 + tstep(1)) * Fsi) .^ 1, ylims)
%                         else
%                             imagesc([-tstep(1), tstep(2)], [1, size(dtc, 1)], dtc .^ 1, ylims)
                        plot(linspace(-tstep(1), tstep(2), size(dtc, 2)), median(dtc .^ 1, 1))
                        xlim([-tstep(1), tstep(2)])
%                         ylim(ylims)
%                             imagesc([-1, 5], [1, size(dtc, 1)], dtc(:, (tstep(1) - 1) * Fsi + 1: (5 + tstep(1)) * Fsi) .^ 1, ylims)
%                         end
                        hold on
                        plot([0, 0], get(gca, 'ylim'), 'r', 'linewidth', 1)
                        hold off
                        colormap('jet')
                    end
                end
            end
            
        case 'whisk_no_whisk_mice_combined_with_stats_air'
            %%% local parameters %%%
            dataswhnwh = datain{1};
            [nc, nw, nwc] = size(dataswhnwh);
            ylims = [0, 1];
            
            %%% main %%%
            figure(gcf)
            clf
            i = 1;
            j = 1;
            dtc = [];
            for k = 1: nwc
                tmp = dataswhnwh{i, j, k}{flag}{1};
                tmp = mean(tmp, 3);
                tmp = tmp - mean(tmp(:, 1: tstep(1) * Fsi), 2);
                dtc = [dtc; tmp];
            end
            
            for ii = 1: size(dtc, 1)
                dtc(ii, :) = normalize(dtc(ii, :));
            end
            
            [~, id] = max(dtc, [], 2);
            [~, id] = sort(id);
            dtc = dtc(id, :);
            imagesc([-1, 2], [1, size(dtc, 1)], dtc(:, (tstep(1) - 1) * Fsi + 1: (2 + tstep(1)) * Fsi), ylims)
            hold on
            plot([0, 0], get(gca, 'ylim'), 'r', 'linewidth', 1)
            hold off
            colormap('parula')
            
            
        case 'average_trace_trial_neuron'
            vfs = datain{1};
            heats = datain{2};
            stp = 1 / 20;
            
            %%% main %%%
            figure(gcf)
            clf

            subplot(1, 2, 1)
            xrg = -stp: stp: (2 + stp);
            hold on
            for i = 1: size(vfs, 2)
                plot(xrg, vfs(39: 81, i), 'linewidth', 2)
            end
            hold off
            xlim([-stp, 2 + stp]);
            legend({'air whisker', 'air no whisker', 'vf whisker + whisk', 'vf whisker + no whisk', 'vf no whisker + whisk', 'vf no whisker + no whisk'})

            subplot(1, 2, 2)
            xrg = -stp: stp: (5 + stp);
            hold on
            for i = 1: size(vfs, 2)
                plot(xrg, heats(39: 141, i), 'linewidth', 2)
            end
            hold off
            xlim([-stp, 5 + stp]);
            legend({'air whisker', 'air no whisker', 'heat whisker + whisk', 'heat whisker + no whisk', 'heat no whisker + whisk', 'heat no whisker + no whisk'})
            
        case 'average_trace_trial_neuron_wipe'
            vfs = datain{1};
            heats = datain{2};
            wps = datain{3};
            tstep = datain{4};
            sk = 1;

            %%% main %%%
            figure(gcf)
            clf
            xrg = linspace(-tstep(1), tstep(2), Fsi * sum(tstep) + 1);
            diss = xrg(2) - xrg(1);

            %%% vf %%%
            subplot(1, 2, 1)
            hold on
            for i = 1: size(vfs, 2)
                tmp = conv(vfs(:, i), gausswin(sk) / sum(gausswin(sk)), 'same');
                plot(xrg, tmp, 'linewidth', 2)
            end
            rg = get(gca, 'ylim');
            dis = diff(rg);
            disuse = dis / 2;
            nn = sum(sum(squeeze(cellfun(@length, wps(2, :, :)))));
            rgn = linspace(rg(2), rg(2) + disuse, nn + 1);
            count = 1;
            ccount = 1;
            for i = 1: size(wps, 2)
                for j = 1: size(wps, 3)
                    for k = 1: length(wps{2, i, j})
                        plot([wps{2, i, j}(k), wps{2, i, j}(k)] / Fsi - tstep(1), [rgn(count), rgn(count + 1)], 'color', clrs(ccount, :), 'linewidth', 2)
                        count = count + 1;
                    end
                    ccount = ccount + 1;
                end
            end
            hold off
            xlim([-diss, 2 + diss])
            ylim([rg(1), rgn(end)])
            legend({'air whisker', 'air no whisker', 'vf whisker + whisk', 'vf whisker + no whisk', 'vf no whisker + whisk', 'vf no whisker + no whisk'})

            %%% heat %%%
            subplot(1, 2, 2)
            hold on
            for i = 1: size(vfs, 2)
                tmp = conv(heats(:, i), gausswin(sk) / sum(gausswin(sk)), 'same');
                plot(xrg, tmp, 'linewidth', 2)
            end
            rg = get(gca, 'ylim');
            dis = diff(rg);
            disuse = dis / 2;
            nn = sum(sum(squeeze(cellfun(@length, wps(1, :, :)))));
            rgn = linspace(rg(2), rg(2) + disuse, nn + 1);
            count = 1;
            ccount = 1;
            for i = 1: size(wps, 2)
                for j = 1: size(wps, 3)
                    for k = 1: length(wps{1, i, j})
                        plot([wps{1, i, j}(k), wps{1, i, j}(k)] / Fsi - tstep(1), [rgn(count), rgn(count + 1)], 'color', clrs(ccount, :), 'linewidth', 2)
                        count = count + 1;
                    end
                    ccount = ccount + 1;
                end
            end
            hold off
            ylim([rg(1), rgn(end)])
            xlim([-diss, 5 + diss])
            legend({'air whisker', 'air no whisker', 'heat whisker + whisk', 'heat whisker + no whisk', 'heat no whisker + whisk', 'heat no whisker + no whisk'}, 'location', 'best')
            
        case 'show_STA_neuron'
            sta = datain{1};
            tstep = [-2, 2];
            crg = [0, 0.00002];
            %%% local parameters %%%
            [nc, nw, nm] = size(sta);
            nwc = length(sta{1, 1, 1});
            
            %%% main %%%
            figure(gcf)
            clf
            mn = zeros(nc, nw, nm, nwc);
            mx = zeros(nc, nw, nm, nwc);
            h = cell(nc, nw, nm, nwc);
            for i = 1: nc
                for j = 1: nw
                    for k = 1: nm
                        for kk = 1: nwc
                            h{i, j, k, kk} = subplot(nc * nw, nm * 3, (k - 1) * 3 + kk + ((i - 1) * nw + j - 1) * nm * 3);
                            tp = sta{i, j, k}{kk};
                            ttmp = squeeze(mean(tp, 3));
                            if isvector(ttmp)
                                ttmp = ttmp(:)';
                            end
                            idx = 1: size(ttmp, 1);
%                             [~, idx] = max(ttmp, [], 2);
%                             [~, idx] = sort(idx);
                            mn(i, j, k, kk) = prctile(ttmp(:), 1);
                            mx(i, j, k, kk) = prctile(ttmp(:), 99);
                            imagesc([tstep(1), tstep(2)], [1, length(idx)], ttmp(idx, :))
                            hold on
                            rgy = get(gca, 'ylim');
                            plot([0, 0], rgy, 'r')
                            hold off
                        end
                    end
                end
            end
            colormap('jet')
            for k = 1: nm
                mnc = min(min(min(squeeze(mn(:, :, k, :)))));
                mxc = max(max(max(squeeze(mx(:, :, k, :)))));
                for i = 1: nc
                    for j = 1: nw
                        for kk = 1: nwc
%                             caxis(h{i, j, k, kk}, [mnc, mxc])
                            caxis(h{i, j, k, kk}, crg)
                        end
                    end
                end
            end
            
        case 'show_STA_group_average'
            vfs = datain{1};
            heats = datain{2};
            if length(datain) < 3
                flag = 1;
            else
                flag = datain{3};
            end
            sk = 1;
            xrg = linspace(-2, 2, size(vfs, 1));
            
            %%% main %%%
            figure(gcf)
            clf

            subplot(1, 2, 1)
            hold on
            for i = 1: size(vfs, 2)
                if flag == 1
%                     plot(xrg, vfs(:, i), 'linewidth', 2)
%                     tmp = smooth(vfs(:, i), sk);
                    tmp = conv(vfs(:, i), gausswin(sk) / sum(gausswin(sk)), 'same');
                    plot(xrg, tmp, 'linewidth', 2)
                else
%                     plot(xrg, vfs(:, i) - vfs(37, i), 'linewidth', 2)
%                     tmp = smooth(vfs(:, i), sk);
                    tmp = conv(vfs(:, i), gausswin(sk) / sum(gausswin(sk)), 'same');
                    plot(xrg, tmp - tmp(40), 'linewidth', 2)
                end
            end
            plot([0, 0], get(gca, 'ylim'), 'r')
            hold off
            legend({'vf whisker + whisk', 'vf whisker + no whisk', 'vf no whisker + whisk', 'vf no whisker + no whisk'}, 'location', 'best')

            subplot(1, 2, 2)
            hold on
            for i = 1: size(vfs, 2)
                if flag == 1
%                     plot(xrg, heats(:, i), 'linewidth', 2)
%                     tmp = smooth(heats(:, i), sk);
                    tmp = conv(heats(:, i), gausswin(sk) / sum(gausswin(sk)), 'same');
                    plot(xrg, tmp, 'linewidth', 2)
                else
                    %                     plot(xrg, heats(:, i) - heats(20, i), 'linewidth', 2)
                    %                     tmp = smooth(heats(:, i), sk);
                    tmp = conv(heats(:, i), gausswin(sk) / sum(gausswin(sk)), 'same');
                    plot(xrg, tmp - tmp(40), 'linewidth', 2)
                end
            end
            plot([0, 0], get(gca, 'ylim'), 'r')
            hold off
            legend({'heat whisker + whisk', 'heat whisker + no whisk', 'heat no whisker + whisk', 'heat no whisker + no whisk'}, 'location', 'best')
            
        case 'show_tracked_neuron_demo'
            strct = datain{1};
            pixh = 200;
            pixw = 320;
            
            a = strct;
            roifn = full(a.spatial_footprints_corrected{1}');
            imga = reshape(max(roifn, [], 2), pixh, pixw);
            imgb = reshape(max(full(a.spatial_footprints_corrected{4}), [], 1), 200, 320);
            
            id = a.cell_to_index_map(:, [1, 4]);
            ida = id(sum(id > 0, 2) == 2);
            tmp = a.centroid_locations_corrected{1}(ida, :);
            roifn = roifn(:, ida);
            x = tmp(:, 1);
            y = tmp(:, 2);
            ids = 1: length(x);
            cntrs = cell(1, length(ids));
            thres = 0.8;
            for i = 1: length(ids)
                tmp = full(reshape(roifn(:, ids(i)), pixh, pixw));
                tmp = imgaussfilt(tmp, 3);
                lvl = max(max(tmp)) * thres;
                cntrs{i} = contour(flipud(tmp), [lvl, lvl]);
                cntrs{i} = [cntrs{i}(:, 2: end - 1), cntrs{i}(:, 2)];
            end
            
            %%% plot %%%
            figure(gcf)
            clf
            imshowpair(imga, imgb)
            hold on
            for i = 1: length(ids)
                plot(cntrs{i}(1, :), pixh - cntrs{i}(2, :), 'r')
                text(x(i), y(i), num2str(i), 'color', [1, 1, 1])
            end
            hold off

    
        case 'show_active_neuron_trial_number'
            vfs = datain{1};
            heats = datain{2};
            tpvm = cellfun(@mean, vfs);
            tphm = cellfun(@mean, heats);
            tpvs = cellfun(@(x) std(x) / (2 * sqrt(length(x))), vfs);
            tphs = cellfun(@(x) std(x) / (2 * sqrt(length(x))), heats);
            
            clf
%             subplot(1,2,1)
            hold on
            for i = 1: length(vfs)
                plot(0.3 * (rand(1, length(vfs{i})) - 0.5) + i, vfs{i}, '.', 'color', clrs(1, :), 'markersize', 6)
            end
            for i = 3: length(heats)
                plot(0.3 * (rand(1, length(heats{i})) - 0.5) + i, heats{i}, '.', 'color', clrs(2, :), 'markersize', 6)
            end
            h1 = errorbar(1: 6, tpvm, tpvs, '-o', 'color', clrs(1, :), 'markersize', 6, 'markerfacecolor', clrs(1, :), 'capsize', 8, 'linewidth', 1);
            h2 = errorbar(3: 6, tphm(3: end), tphs(3: end), '-s', 'color', clrs(2, :), 'markersize', 6, 'markerfacecolor', clrs(2, :), 'capsize', 8, 'linewidth', 1);
            hold off
            xticks(1: 6)
            xticklabels({'air w', 'air n', 'w w', 'w nw', 'n w', 'n nw'})
            legend([h1, h2], {'vf', 'heat'}, 'location', 'best')
            ylim([0, 1])
            
%             subplot(1,2,2)
%             hold on
%             plot(vfs(1,:) ./ vfs(2,:))
%             plot(heats(1,:) ./ heats(2,:))
%             hold off
%             xticklabels({'air w', 'air n', 'w w', 'w nw', 'n w', 'n nw'})
%             legend({'vf', 'heat'}, 'location', 'best')
            
        case 'show_all_signals'
            sigs = datain{1};
            bidx = datain{2};
            sigs = 2 * sigs([5:9,12,13,23,41,42,49,56,57,62,64,67:69,71,72,74,78,82,84], 1201: 5: 5200);
            bidx = bidx(1201: 5: 5200);
% %             tmp = std(sigs, [], 2);
% %             [~, iduse] = sort(tmp, 'descend');
% %             sigs = sigs(iduse, :);
% %             for i = 1: size(sigs, 1)
% %                 sigs(i, :) = 0.1 * normalize(sigs(i, :));
% %             end
            
            b = sigs;
            scl = 0.1;
            figure(gcf)
            clf
            hold on
            plot((b + scl * (1:size(b,1))')','k');
            axis tight;
            rg = get(gca, 'ylim');
            [l, n] = bwlabeln(bidx);
            for ii = 1: n
                tmp = find(l == ii);
%                 plot(tmp, repmat(rg(2) + scl, 1, length(tmp)), 'g', 'linewidth', 2)
                plot([tmp(1), tmp(1)], rg, 'g')
            end
            ylim([0, 0.1 * size(sigs, 1) + 0.2])
            
        case 'sct_neuron_mean_std_event_vis'
            stats = datain{1};
            distr = datain{2};
            srgy = datain{3};
            drgy = datain{4};
            flag = datain{5};
            [nc, nw, nm] = size(stats);
            
            %%% plot 3 conds, all mice in one plot but diff colors %%%
            figure(gcf)
            clf
            for i = 1: nc
                subplot(3, 3, i)
                hold on
                errorbar(stats{i, 1}{flag}(:, 3), stats{i, 1}{flag}(:, 1), stats{i, 1}{flag}(:, 2) / 2, 'd', 'color', clrs(2 * (i - 1) + 1, :), 'markerfacecolor', clrs(2 * (i - 1) + 1, :), 'markersize', 2)
                errorbar(stats{i, 2}{flag}(:, 3), stats{i, 2}{flag}(:, 1), stats{i, 2}{flag}(:, 2) / 2, 's', 'color', clrs(2 * i, :), 'markerfacecolor', clrs(2 * i, :), 'markersize', 2)
                hold off
                ylim(srgy(flag, :))
                title(opt.conds{i})
            end
            
            %%% plot event frequency dist %%%
            for i = 1: nc
                subplot(3, 3, 3 + i)
                hold on
                tmp = distr{i, 1}{flag}(:, 1);
                plot(linspace(0, 1, length(tmp)), tmp, 'linewidth', 2, 'color', clrs(2 * (i - 1) + 1, :))
                tmp = distr{i, 2}{flag}(:, 1);
                plot(linspace(0, 1, length(tmp)), tmp, 'linewidth', 2, 'color', clrs(2 * i, :))
                hold off
                xlim([0, 1])
                ylim(drgy{flag}(1, :))
            end
            
            %%% plot mean amplitude dist %%%
            for i = 1: nc
                subplot(3, 3, 6 + i)
                hold on
                tmp = distr{i, 1}{flag}(:, 2);
                plot(linspace(srgy(flag, 1), srgy(flag, 2), length(tmp)), tmp, 'linewidth', 2, 'color', clrs(2 * (i - 1) + 1, :))
%                 plot(linspace(srgy(flag, 1), srgy(flag, 2), length(tmp)), cumsum(tmp), 'linewidth', 2, 'color', clrs(2 * (i - 1) + 1, :))
                tmp = distr{i, 2}{flag}(:, 2);
                plot(linspace(srgy(flag, 1), srgy(flag, 2), length(tmp)), tmp, 'linewidth', 2, 'color', clrs(2 * i, :))
%                 plot(linspace(srgy(flag, 1), srgy(flag, 2), length(tmp)), cumsum(tmp), 'linewidth', 2, 'color', clrs(2 * i, :))
                hold off
                xlim(srgy(flag, :) / 3)
                ylim(drgy{flag}(2, :))
            end
            
        case 'sct_neuron_mean_std_event_whnwh_vis'
            stats = datain{1};
            distr = datain{2};
            srgy = datain{3};
            drgy = datain{4};
            flag = datain{5};
            [nc, nw, nm] = size(stats);
            
            %%% plot 3 conds, all mice in one plot but diff colors %%%
            figure(gcf)
            clf
            for i = 1: nc
                subplot(3, 6, 2 * (i - 1) + 1)
                hold on
                errorbar(stats{i, 1}{flag}{1}(:, 3), stats{i, 1}{flag}{1}(:, 1), stats{i, 1}{flag}{1}(:, 2) / 2, 'd', 'color', clrs(2 * (i - 1) + 1, :), 'markerfacecolor', clrs(2 * (i - 1) + 1, :), 'markersize', 2)
                errorbar(stats{i, 2}{flag}{1}(:, 3), stats{i, 2}{flag}{1}(:, 1), stats{i, 2}{flag}{1}(:, 2) / 2, 's', 'color', clrs(2 * i, :), 'markerfacecolor', clrs(2 * i, :), 'markersize', 2)
                hold off
                ylim(srgy(flag, :))
                title([opt.conds{i}, ', whisk'])
                legend({'Whisker', 'No Whisker'})
            end
            for i = 1: nc
                subplot(3, 6, 2 * i)
                hold on
                errorbar(stats{i, 1}{flag}{2}(:, 3), stats{i, 1}{flag}{2}(:, 1), stats{i, 1}{flag}{2}(:, 2) / 2, 'd', 'color', clrs(2 * (i - 1) + 1, :), 'markerfacecolor', clrs(2 * (i - 1) + 1, :), 'markersize', 2)
                errorbar(stats{i, 2}{flag}{2}(:, 3), stats{i, 2}{flag}{2}(:, 1), stats{i, 2}{flag}{2}(:, 2) / 2, 's', 'color', clrs(2 * i, :), 'markerfacecolor', clrs(2 * i, :), 'markersize', 2)
                hold off
                ylim(srgy(flag, :))
                title([opt.conds{i}, ', no whisk'])
            end
            
            %%% plot event frequency dist %%%
            for i = 1: nc
                subplot(3, 6, 6 + 2 * (i - 1) + 1)
                hold on
                tmp = distr{i, 1}{flag}{1}(:, 1);
                plot(linspace(0, 1, length(tmp)), tmp, 'linewidth', 2, 'color', clrs(2 * (i - 1) + 1, :))
                tmp = distr{i, 2}{flag}{1}(:, 1);
                plot(linspace(0, 1, length(tmp)), tmp, 'linewidth', 2, 'color', clrs(2 * i, :))
                hold off
                xlim([0, 1])
                ylim(drgy{flag}(1, :))
            end
            for i = 1: nc
                subplot(3, 6, 6 + 2 * i)
                hold on
                tmp = distr{i, 1}{flag}{2}(:, 1);
                plot(linspace(0, 1, length(tmp)), tmp, 'linewidth', 2, 'color', clrs(2 * (i - 1) + 1, :))
                tmp = distr{i, 2}{flag}{2}(:, 1);
                plot(linspace(0, 1, length(tmp)), tmp, 'linewidth', 2, 'color', clrs(2 * i, :))
                hold off
                xlim([0, 1])
                ylim(drgy{flag}(1, :))
            end
            
            %%% plot mean amplitude dist %%%
            for i = 1: nc
                subplot(3, 6, 12 + 2 * (i - 1) + 1)
                hold on
                tmp = distr{i, 1}{flag}{1}(:, 2);
                plot(linspace(srgy(flag, 1), srgy(flag, 2), length(tmp)), tmp, 'linewidth', 2, 'color', clrs(2 * (i - 1) + 1, :))
                tmp = distr{i, 2}{flag}{1}(:, 2);
                plot(linspace(srgy(flag, 1), srgy(flag, 2), length(tmp)), tmp, 'linewidth', 2, 'color', clrs(2 * i, :))
                hold off
                xlim(srgy(flag, :) / 3)
                ylim(drgy{flag}(2, :))
            end
            for i = 1: nc
                subplot(3, 6, 12 + 2 * i)
                hold on
                tmp = distr{i, 1}{flag}{2}(:, 2);
                plot(linspace(srgy(flag, 1), srgy(flag, 2), length(tmp)), tmp, 'linewidth', 2, 'color', clrs(2 * (i - 1) + 1, :))
                tmp = distr{i, 2}{flag}{2}(:, 2);
                plot(linspace(srgy(flag, 1), srgy(flag, 2), length(tmp)), tmp, 'linewidth', 2, 'color', clrs(2 * i, :))
                hold off
                xlim(srgy(flag, :) / 3)
                ylim(drgy{flag}(2, :))
            end
                   
        case 'sct_neuron_mean_std_event_whnwh_vis_only_full'
            stats = datain{1};
            distr = datain{2};
            srgy = datain{3};
            drgy = datain{4};
            flag = datain{5};
            [nc, nw, nm] = size(stats);
            
            %%% plot 3 conds, all mice in one plot but diff colors %%%
            figure(gcf)
            clf
            for i = 1: nc
                subplot(3, 3, i)
                hold on
                errorbar(stats{i, 1}{flag}{1}(:, 3), stats{i, 1}{flag}{1}(:, 1), stats{i, 1}{flag}{1}(:, 2) / 2, 'd', 'color', clrs(2 * (i - 1) + 1, :), 'markerfacecolor', clrs(2 * (i - 1) + 1, :), 'markersize', 2)
                errorbar(stats{i, 1}{flag}{2}(:, 3), stats{i, 1}{flag}{2}(:, 1), stats{i, 1}{flag}{2}(:, 2) / 2, 's', 'color', clrs(2 * i, :), 'markerfacecolor', clrs(2 * i, :), 'markersize', 2)
                hold off
                ylim(srgy(flag, :))
                title([opt.conds{i}, ', whisk'])
                legend({'Whisk', 'Non-whisk'})
            end
            
            %%% plot event frequency dist %%%
            for i = 1: nc
                subplot(3, 3, 3 + i)
                hold on
                tmp = distr{i, 1}{flag}{1}(:, 1);
                plot(linspace(0, 1, length(tmp)), tmp, 'linewidth', 2, 'color', clrs(2 * (i - 1) + 1, :))
                tmp = distr{i, 1}{flag}{2}(:, 1);
                plot(linspace(0, 1, length(tmp)), tmp, 'linewidth', 2, 'color', clrs(2 * i, :))
                hold off
                xlim([0, 1])
                ylim(drgy{flag}(1, :))
            end
            
            %%% plot mean amplitude dist %%%
            for i = 1: nc
                subplot(3, 3, 6 + i)
                hold on
                tmp = distr{i, 1}{flag}{1}(:, 2);
                plot(linspace(srgy(flag, 1), srgy(flag, 2), length(tmp)), tmp, 'linewidth', 2, 'color', clrs(2 * (i - 1) + 1, :))
                tmp = distr{i, 1}{flag}{2}(:, 2);
                plot(linspace(srgy(flag, 1), srgy(flag, 2), length(tmp)), tmp, 'linewidth', 2, 'color', clrs(2 * i, :))
                hold off
                xlim(srgy(flag, :) / 3)
                ylim(drgy{flag}(2, :))
            end
            
        case 'sct_neuron_mean_std_event_whnwh_vis_only_full_all_trials'
            dataextswhnwhaf = datain{1};
            flag = datain{2};
            [nc, nw, nm] = size(dataextswhnwhaf);
            c = cell(3, 2);
            for i = 1: nc
                for j = 1: 1
                    for k = 1: nm
                        for jj = 1: 2
                            tmp = dataextswhnwhaf{i, j, k}{flag}{jj};
                            tmp1 = squeeze(max(tmp, [], 2));
%                             idx = tmp1 > 2 * std(tmp(:));
%                             tmp2 = tmp1(idx);
%                             c{i, jj} = [c{i, jj}; tmp2(:)];
                            c{i, jj} = [c{i, jj}; tmp1(:)];
                        end
                    end
                end
            end
            
            edges = linspace(min(cell2mat(c(:))), max(cell2mat(c(:))), 101);
            ctrs = (edges(1: end - 1) + edges(2: end)) / 2;
            d = cell(nc, 2);
            for i = 1: nc
                for j = 1: 2
                    tmp = c{i, j};
%                     d{i, j} = ksdensity(tmp, ctrs, 'bandwidth', (ctrs(end) - ctrs(1)) / (5 * length(ctrs)));
                    d{i, j} = histcounts(tmp, edges);
                end
            end
            
            figure(gcf)
            clf
            hold on
            plot(ctrs, d{2,1} / sum(d{2, 1}))
            plot(ctrs, d{2,2} / sum(d{2, 2}))
            plot(ctrs, d{3,1} / sum(d{3, 1}))
            plot(ctrs, d{3,2} / sum(d{3, 2}))
%             plot(ctrs, cumsum(d{2, 1}) / sum(d{2, 1}))
%             plot(ctrs, cumsum(d{2, 2}) / sum(d{2, 2}))
%             plot(ctrs, cumsum(d{3, 1}) / sum(d{3, 1}))
%             plot(ctrs, cumsum(d{3, 2}) / sum(d{3, 2}))
            
        case 'sct_neuron_mean_std_event_combined_vis'
            stats = datain{1};
            distr = datain{2};
            srgy = datain{3};
            drgy = datain{4};
            flag = datain{5};
            [nc, nw, nm] = size(stats);
            
            %%% plot event frequency dist %%%
            figure(gcf)
            clf
            subplot(2, 1, 1)
            hold on
            for i = 1: nc
                tmp = distr{i, 1}{flag}(:, 1);
                plot(linspace(0, 1, length(tmp)), tmp, 'linewidth', 2, 'color', clrs(2 * (i - 1) + 1, :))
                tmp = distr{i, 2}{flag}(:, 1);
                plot(linspace(0, 1, length(tmp)), tmp, 'linewidth', 2, 'color', clrs(2 * i, :))
                xlim([0, 1])
                ylim(drgy{flag}(1, :))
            end
            hold off
            lgds1 = repmat(opt.conds, 2, 1);
            lgds2 = repmat(opt.whisks', 1, 3);
            lgds = cellfun(@(x, y) [x, ', ', y], lgds1, lgds2, 'uniformoutput', false);
            legend(lgds(:))
            
            %%% plot mean amplitude dist %%%
            subplot(2, 1, 2)
            for i = 1: nc
                hold on
                tmp = distr{i, 1}{flag}(:, 2);
                plot(linspace(srgy(flag, 1), srgy(flag, 2), length(tmp)), tmp, 'linewidth', 2, 'color', clrs(2 * (i - 1) + 1, :))
%                 plot(linspace(srgy(flag, 1), srgy(flag, 2), length(tmp)), cumsum(tmp), 'linewidth', 2, 'color', clrs(2 * (i - 1) + 1, :))
                tmp = distr{i, 2}{flag}(:, 2);
                plot(linspace(srgy(flag, 1), srgy(flag, 2), length(tmp)), tmp, 'linewidth', 2, 'color', clrs(2 * i, :))
%                 plot(linspace(srgy(flag, 1), srgy(flag, 2), length(tmp)), cumsum(tmp), 'linewidth', 2, 'color', clrs(2 * i, :))
                hold off
                xlim(srgy(flag, :) / 3)
                ylim(drgy{flag}(2, :))
            end
            
        case 'sct_neuron_mean_std_event_whnwh_cumsum_vis'
            stats = datain{1};
            distr = datain{2};
            srgy = datain{3};
            drgy = datain{4};
            flag = datain{5};
            [nc, nw] = size(stats);
            rgy = srgy(flag, :);
            edges = linspace(rgy(1), rgy(2), 41);
            ctrs = (edges(1: end - 1) + edges(2: end)) / 2;
            
            figure(gcf)
            clf
            hold on
            for i = 2: nc
                for j = 1: nw
                    tmp = stats{i, j}{flag}{1}(:, 1);
                    t = histcounts(tmp, edges);
                    if j == 1
                        plot(ctrs, cumsum(t), 'linewidth', 2)
                    else
                        plot(ctrs, cumsum(t))
                    end
                    
                    tmp = stats{i, j}{flag}{2}(:, 1);
                    t = histcounts(tmp, edges);
                    plot(ctrs, cumsum(t))
                end
            end
            legend({'h.w.w.', 'h.w.nw.', 'h.nw.w.', 'h.nw.nw', 'v.w.w.', 'v.w.nw.', 'v.nw.w.', 'v.nw.nw'})
            
        case 'sct_neuron_mean_std_event_whnwh_high_vis'
            stats = datain{1};
            distr = datain{2};
            srgy = datain{3};
            drgy = datain{4};
            flag = datain{5};
            [nc, nw] = size(stats);
            rgy = srgy(flag, :);
            edges = linspace(rgy(1), rgy(2), 41);
            ctrs = (edges(1: end - 1) + edges(2: end)) / 2;
            thres = ctrs(4);
            
            dt = [];
            for i = 2: nc
                for j = 1: nw
                    for k = 1: 2
                        tmp = stats{i, j}{flag}{k}(:, 1);
                        dt = [dt, sum(tmp > thres)];
                    end
                end
            end
            
            figure(gcf)
            clf
            hold on
            dtfw = dt([1, 5]);
            dtr = dt(setdiff(1: 8, [1, 5]));
            bar(1: length(dtfw), dtfw)
            bar((1: length(dtr)) + length(dtfw), dtr)
            
            
            
            
        case 'sct_neuron_mean_std_event_whnwh_combined_vis'
            stats = datain{1};
            distr = datain{2};
            srgy = datain{3};
            drgy = datain{4};
            flag = datain{5};
            [nc, nw, nm] = size(stats);
            
            %%% plot event frequency dist %%%
            figure(gcf)
            clf
            subplot(2, 1, 1)
            hold on
            for i = 2: nc
                tmp = distr{i, 1}{flag}{1}(:, 1);
                plot(linspace(0, 1, length(tmp)), tmp, 'linewidth', 2, 'color', clrs(2 * (i - 2) + 1, :))
                tmp = distr{i, 2}{flag}{1}(:, 1);
                plot(linspace(0, 1, length(tmp)), tmp, 'linewidth', 2, 'color', clrs(2 * (i - 1), :))
                xlim([0, 1])
                ylim(drgy{flag}(1, :))
            end
            for i = 2: nc
                tmp = distr{i, 1}{flag}{2}(:, 1);
                plot(linspace(0, 1, length(tmp)), tmp, 'linewidth', 2, 'color', clrs(4 + 2 * (i - 2) + 1, :))
                tmp = distr{i, 2}{flag}{2}(:, 1);
                plot(linspace(0, 1, length(tmp)), tmp, 'linewidth', 2, 'color', clrs(4 + 2 * (i - 1), :))
                xlim([0, 1])
                ylim(drgy{flag}(1, :))
            end
            hold off
            lgds1 = repmat(opt.conds(2: end), 2, 2);
            lgds2 = repmat(opt.whisks', 1, 4);
            lgds3 = repmat({'whisk', 'no whisk'}, 2, 2);
            lgds = cellfun(@(x, y, z) [x, ', ', y, ', ', z], lgds1, lgds2, lgds3, 'uniformoutput', false);
            legend(lgds(:))
            
            %%% plot mean amplitude dist %%%
            subplot(2, 1 ,2)
            hold on
            for i = 2: nc
                tmp = distr{i, 1}{flag}{1}(:, 2);
                plot(linspace(srgy(flag, 1), srgy(flag, 2), length(tmp)), tmp, 'linewidth', 2, 'color', clrs(2 * (i - 1) + 1, :))
                tmp = distr{i, 2}{flag}{1}(:, 2);
                plot(linspace(srgy(flag, 1), srgy(flag, 2), length(tmp)), tmp, 'linewidth', 2, 'color', clrs(2 * i, :))
                xlim(srgy(flag, :) / 3)
                ylim(drgy{flag}(2, :))
            end
            for i = 2: nc
                tmp = distr{i, 1}{flag}{2}(:, 2);
                plot(linspace(srgy(flag, 1), srgy(flag, 2), length(tmp)), tmp, 'linewidth', 2, 'color', clrs(2 * (i - 1) + 1, :))
                tmp = distr{i, 2}{flag}{2}(:, 2);
                plot(linspace(srgy(flag, 1), srgy(flag, 2), length(tmp)), tmp, 'linewidth', 2, 'color', clrs(2 * i, :))
                xlim(srgy(flag, :) / 3)
                ylim(drgy{flag}(2, :))
            end
            hold off

        case 'sct_same_neuron_compare_vis'
            data = datain{1};
            flag = datain{2};
            mn = cell2mat(cellfun(@(x) mean(log10(x{flag}(~isnan(x{flag}(:, 1)) & sum(abs(x{flag}) == Inf | x{flag} == 0, 2) == 0, :)), 1), data, 'uniformoutput', false)');
            ste = sqrt(std(mn, [], 1) / size(mn, 1));
            tp = cell2mat(cellfun(@(x) x{flag}, data, 'uniformoutput', false)');
            tp = log10(tp(~isnan(tp(:, 1)) & sum(abs(tp) == Inf | tp == 0, 2) == 0, :));
            mn1 = mean(tp, 1);
            ste1 = sqrt(std(tp, [], 1) / size(tp, 1));
            
            figure(gcf)
            clf
            hold on
            errorbar(mean(mn, 1), ste)
            errorbar(mn1, ste1)
            xlim([0, size(mn, 2) + 1])
            xticks(1: size(mn, 2))
            lgds1 = repmat(opt.conds, 2, 1)';
            lgds2 = repmat(opt.whisks', 1, 3)';
            lgds = cellfun(@(x, y) [x, ', ', y], lgds1, lgds2, 'uniformoutput', false);
            xticklabels(lgds)
            legend({'mouse average', 'all neuron average'})
            
        case 'sct_same_neuron_compare_ratio_stats_vis'
            data = datain{1};
            idx = datain{2};
            flag = datain{3};
            name1 = {'heat full', 'vf full', 'air no', 'heat no', 'vf no'};
            name2 = {'inc', 'dec'};
            names = cellfun(@(x, y) [x, ', ', y], repmat(name1, 2, 1), repmat(name2', 1, 5), 'uniformoutput', false);
            
            figure(gcf)
            clf
            subplot(2, 2, 1)
            pie(data{1, flag})
            
            subplot(2, 2, 2)
            tmp = data{1, flag};
            idt = idx{1, flag};
            idref = [1, 0];
            tp1 = zeros(2, 5);
            for i = 1: 2
                for j = 1: 5
                    tp1(i, j) = sum(tmp(idt(:, j) == idref(i)));
                end
            end
            pie(tp1(:), names)
            
            subplot(2, 2, 3)
            pie(data{2, flag})
            
            subplot(2, 2, 4)
            tmp = data{2, flag};
            idt = idx{2, flag};
            idref = [1, 0];
            tp2 = zeros(2, 5);
            for i = 1: 2
                for j = 1: 5
                    tp2(i, j) = sum(tmp(idt(:, j) == idref(i)));
                end
            end
            pie(tp2(:), names)
        
        case 'sct_condition_distribution_vis'
            sct_air = datain{1};
            sct_an = datain{2};
            
            tp1 = 100 * sct_air ./ sct_air(1, :);
            tp2 = 100 * sct_air ./ sct_an;
            tpp1 = 100 * sum(sct_air, 2) ./ sum(sct_air(1, :));
            tpp2 = 100 * sum(sct_air, 2) ./ sum(sct_an, 2);
            xrg1 = 0.8: 5.8;
            xrg2 = 1.2: 6.2;
            
            figure(gcf)
            clf
            hold on
            p1 = plot(xrg1, tp1, 'color', [0.5, 0.5, 0.5], 'linewidth', 5);
            p2 = plot(xrg2, tp2, 'color', [0.4, 0.7, 0.7], 'linewidth', 5);
            for j = 1: length(p1)
                p1(j).Color(4) = 0.4;
                p2(j).Color(4) = 0.4;
            end
            plot(xrg1, tpp1, 'color', [0.5, 0.5, 0.5], 'linewidth', 10)
            plot(xrg2, tpp2, 'color', [0.4, 0.7, 0.7], 'linewidth', 10)
%             plot(xrg1, tpp1, '.k', 'linewidth', 10)
%             plot(xrg2, tpp2, '.k', 'linewidth', 10)
%             b = bar([tpp1, tpp2], 'facecolor', 'flat');
            b1 = bar(xrg1, tpp1, 0.3, 'facecolor', 'flat');
            b2 = bar(xrg2, tpp2, 0.3, 'facecolor', 'flat');
            b1.CData = repmat([0.5, 0.5, 0.5], size(b1.CData, 1), 1);
            b2.CData = repmat([0.4, 0.7, 0.7], size(b2.CData, 1), 1);
            errorbar(xrg1, tpp1, std(tp1, [], 2) / 2, '.k', 'markersize', 15, 'linewidth', 2, 'capsize', 8)
            alpha(0.8)
            errorbar(xrg2, tpp2, std(tp2, [], 2) / 2, '.', 'color', [0.3, 0.6, 0.6], 'markersize', 15, 'linewidth', 2, 'capsize', 8)
            alpha(0.8)
            legend([b1, b2], {'% active/air full', '% active/each condition'})
%             xlim([0, 1])
            ylim([0, 120])
            yticks(0: 20: 100)
            xticks(1: 6)
            xticklabels({'air f', 'heat f', 'vf f', 'air no', 'heat no', 'vf no'})
            ylabel('Percentage')
            
            
            
        % --------------------------------------------- %
        
        case 'sct_condition_distribution_overlap_full_vis'
            dt = datain{1};
            conds = 1: 7;
            nm = length(dt);
            
            tp = zeros(1, length(conds));
            for i = 1: nm
                for j = conds
                    tmp = dt{i};
                    tmp = tmp(:, 1: 2: 5);
                    ref = de2bi(conds(j), 3);
                    tp(j) = tp(j) + sum(ismember(tmp, ref, 'rows'));
                end
            end
            
            figure(gcf)
            clf
            lbl = {'a', 'h', 'v', 'ah', 'av', 'hv', 'ahv'};
            tp = tp([4, 2, 1, 6, 5, 3, 7]);
            plot(conds, tp);
            for i = conds
                text(i, tp(i) + 1, num2str(tp(i)))
            end
            xticklabels(lbl)
            
        case 'sct_condition_distribution_overlap_no_vis'
            dt = datain{1};
            conds = 1: 7;
            nm = length(dt);
            
            tp = zeros(1, length(conds));
            for i = 1: nm
                for j = conds
                    tmp = dt{i};
                    tmp = tmp(:, 2: 2: 6);
                    ref = de2bi(conds(j), 3);
                    tp(j) = tp(j) + sum(ismember(tmp, ref, 'rows'));
                end
            end
            
            figure(gcf)
            clf
            lbl = {'a', 'h', 'v', 'ah', 'av', 'hv', 'ahv'};
            tp = tp([4, 2, 1, 6, 5, 3, 7]);
            plot(conds, tp);
            for i = conds
                text(i, tp(i) + 1, num2str(tp(i)))
            end
            xticklabels(lbl)
            
        case 'tca_3d_vis'
            data = datain{1};
            [nc, nw, nm] = size(data);
            mn = min(cell2mat(cellfun(@(x) min(x, [], 1), data(:), 'uniformoutput', false)), [], 1);
            mx = max(cell2mat(cellfun(@(x) max(x, [], 1), data(:), 'uniformoutput', false)), [], 1);
            wdw = [2 * Fsi + 1, 4 * Fsi; 2 * Fsi, 7 * Fsi; 2 * Fsi, 3 * Fsi];
            
            figure(gcf)
            clf
            for i = 1: nc
                for j = 1: nw
                    subplot(nw, nc, nc * (j - 1) + i)
                    hold on
                    for k = 1: nm
                        tmp = data{i, j, k};
%                         plot3(tmp(wdw(i, 1): wdw(i, 2), 1), tmp(wdw(i, 1): wdw(i, 2), 2), tmp(wdw(i, 1): wdw(i, 2), 3), 'color', clrs(1, :))
%                         plot(tmp(wdw(i, 1): wdw(i, 2), 1), tmp(wdw(i, 1): wdw(i, 2), 2), 'color', clrs(1, :))
%                         plot(tmp(wdw(i, 1): wdw(i, 2), 1), tmp(wdw(i, 1): wdw(i, 2), 2))
                        plot(tmp(wdw(i, 1): wdw(i, 2), 1))
                    end
                    ylim([mn(1), mx(1)])
%                     xlim([mn(1), mx(1)])
%                     ylim([mn(2), mx(2)])
%                     zlim([mn(3), mx(3)])
                    hold off
%                     view(45, 15)
                end
            end
         
        case 'tca_behav_vis'
            tp = datain{1};
            rgy = datain{2};
            [nc, nw] = size(tp);
            
            count = 1;
%             cols = length(tp{1, 1});
            cols = 4;
            xrg = (-tstep(1) * Fsi: tstep(2) * Fsi) / Fsi;
            skp = 4;
            cnames = {'whnwp', 'nwhnwp', 'whwp', 'nwhwp', 'wh', 'nwh', 'wp', 'nwp'};
            rnames = {'heat w', 'heat nw', 'vf w', 'vf nw'};
            
            figure(gcf)
            clf
            for i = 1: nc
                for j = 1: nw
                    tmp = tp{i, j};
                    for k = 1: cols
                        subplot(nc * nw, cols, (count - 1) * cols + k)
%                         tpp = squeeze(tmp(:, :, k));
%                         [~, idt] = max(tpp, [], 2);
%                         [~, idt] = sort(idt);
%                         imagesc(xrg, 1: size(tmp, 1), tpp(idt, :), [0, 0.2])
                        hold on
                        plot(xrg, squeeze(tmp{k})', 'color', clrs(1, :))
                        plot(xrg, squeeze(median(tmp{k}, 1))', 'color', clrs(1, :), 'linewidth', 2)
                        errorbar(xrg(1: skp: end), squeeze(median(tmp{k}(:, 1: skp: end), 1)), squeeze(std(tmp{k}(:, 1: skp: end), [], 1)) / (2 * sqrt(size(tmp{k}, 1))));
%                         plot(xrg, squeeze(mean(tmp{k}, 1))', 'color', clrs(1, :), 'linewidth', 2)
%                         errorbar(xrg(1: skp: end), squeeze(mean(tmp{k}(:, 1: skp: end), 1)), squeeze(std(tmp{k}(:, 1: skp: end), [], 1)) / (2 * sqrt(size(tmp{k}, 1))));
                        xlim([-tstep(1), tstep(2)])
                        ylim(rgy)
                        hold off
                        
                        if k == 1
                            ylabel(rnames{count})
                        end
                        
                        if count == 1
                            title(cnames{k})
                        end
                    end
                    count = count + 1;
                end
            end
            
        case 'tca_behav_vis_heatmap'
            tp = datain{1};
            [nc, nw] = size(tp);
            
            count = 1;
%             cols = length(tp{1, 1});
            cols = 4;
            xrg = (-tstep(1) * Fsi: tstep(2) * Fsi) / Fsi;
            skp = 4;
            cnames = {'whnwp', 'nwhnwp', 'whwp', 'nwhwp', 'wh', 'nwh', 'wp', 'nwp'};
            rnames = {'heat w', 'heat nw', 'vf w', 'vf nw'};
            
            figure(gcf)
            clf
            colormap('jet')
            for i = 1: nc
                for j = 1: nw
                    tmp = tp{i, j};
                    for k = 1: cols
                        subplot(nc * nw, cols, (count - 1) * cols + k)
                        tpp = squeeze(tmp{k});
                        [~, idt] = max(tpp, [], 2);
                        [~, idt] = sort(idt);
                        for ii = 1: size(tpp, 1)
                            tpp(ii, :) = normalize(tpp(ii, :));
                        end
                        imagesc(xrg, 1: size(tpp, 1), imgaussfilt(tpp(idt, :), 1) .^ 1, [0, 1])
%                         imagesc(xrg, 1: size(tpp, 1), tpp(idt, :) .^ 1, [0, 3])
%                         hold on
%                         plot(xrg, squeeze(tmp{k})', 'color', clrs(1, :))
%                         plot(xrg, squeeze(median(tmp{k}, 1))', 'color', clrs(1, :), 'linewidth', 2)
%                         errorbar(xrg(1: skp: end), squeeze(median(tmp{k}(:, 1: skp: end), 1)), squeeze(std(tmp{k}(:, 1: skp: end), [], 1)) / (2 * sqrt(size(tmp, 1))));
%                         xlim([-tstep(1), tstep(2)])
%                         ylim([0, 3])
%                         hold off
                        
                        if k == 1
                            ylabel(rnames{count})
                        end
                        
                        if count == 1
                            title(cnames{k})
                        end
                    end
                    count = count + 1;
                end
            end

            
        case 'tca_behav_vis_subcluster_temp'
            tp = datain{1};
            if length(datain) < 2
                rg = [0, 1];
            else
                rg = datain{2};
            end
            [nc, nw] = size(tp);
            
            count = 1;
%             cols = length(tp{1, 1});
            cols = 4;
            xrg = (-tstep(1) * Fsi: tstep(2) * Fsi) / Fsi;
            skp = 4;
            cnames = {'whnwp', 'nwhnwp', 'whwp', 'nwhwp', 'wh', 'nwh', 'wp', 'nwp'};
            rnames = {'heat w', 'heat nw', 'vf w', 'vf nw'};
            lnames = {'1st [0-1s]', '2nd [1-2s]', '3rd [2-3s]', '4th [3-5s]'};
            
            figure(gcf)
            clf
            for i = 1: nc
                for j = 1: nw
                    tmp = tp{i, j};
                    for k = 1: cols
                        subplot(nc * nw, cols, (count - 1) * cols + k)
                        p = [];
%                         tpp = squeeze(tmp(:, :, k));
%                         [~, idt] = max(tpp, [], 2);
%                         [~, idt] = sort(idt);
%                         imagesc(xrg, 1: size(tmp, 1), tpp(idt, :), [0, 0.2])
                        hold on
                        [nn, nt, nrp] = size(tmp{k});
                        for ii = 1: nn
                            t1 = tmp{k}(ii, :, 1);
                            t2 = tmp{k}(ii, :, 2);
                            x = [xrg, fliplr(xrg)];
                            y = [t1 + t2, fliplr(t1 - t2)];
                            patch(x, 10 * (y + 0.1 * ii), clrs(mod(ii, 8), :), 'edgecolor', 'none')
                            alpha(0.3)
                            p(ii) = plot(xrg, t1, 'color', clrs(mod(ii - 1, 8) + 1, :), 'linewidth', 1);
%                             p(ii) = plot(xrg, 10 * (t1 + 0.1 * ii), 'color', clrs(mod(ii - 1, 8) + 1, :), 'linewidth', 1);
                            
                        end
%                         t1 = squeeze(mean(tmp{k}, 1));
%                         plot(xrg, t1, 'k', 'linewidth', 4)
                        xlim([-tstep(1), tstep(2)])
                        ylim(rg)
                        hold off
                        
                        if k == 1
                            ylabel(rnames{count})
                        end
                        
                        if count == 1
                            title(cnames{k})
                        end
                        
                        if count == 1 && k == 1
                            legend(p, lnames)
                        end
                    end
                    count = count + 1;
                end
            end
 
        case 'tca_behav_vis_subcluster_dist'
            tp = datain{1};
            [nc, nw] = size(tp);
            
            count = 1;
%             cols = length(tp{1, 1});
            cols = 4;
%             xrg = (-tstep(1) * Fsi: tstep(2) * Fsi) / Fsi;
            skp = 4;
            cnames = {'whnwp', 'nwhnwp', 'whwp', 'nwhwp', 'wh', 'nwh', 'wp', 'nwp'};
            rnames = {'heat w', 'heat nw', 'vf w', 'vf nw'};
            
            figure(gcf)
            clf
            for i = 1: nc
                if i == 1
                    xrg = linspace(-2, 10, size(tp{1, 1}, 2) + 1);
                    xrg = (xrg(1: end - 1) + xrg(2: end)) / 2;
                else
                    xrg = linspace(-0.5, 1, size(tp{1, 1}, 2) + 1);
                    xrg = (xrg(1: end - 1) + xrg(2: end)) / 2;
                end
                for j = 1: nw
                    tmp = tp{i, j};
                    for k = 1: cols
                        subplot(nc * nw, cols, (count - 1) * cols + k)
%                         tpp = squeeze(tmp(:, :, k));
%                         [~, idt] = max(tpp, [], 2);
%                         [~, idt] = sort(idt);
%                         imagesc(xrg, 1: size(tmp, 1), tpp(idt, :), [0, 0.2])
                        hold on
%                         [nn, ~, nrp] = size(tmp);
                        plot(xrg, tmp(k, :))
                        if i == 1
                            xlim([-2, 10])
                        else
                            xlim([-0.5, 1])
                        end
                        ylim([0, 0.1])
                        hold off
                        
                        if k == 1
                            ylabel(rnames{count})
                        end
                        
                        if count == 1
                            title(cnames{k})
                        end
                    end
                    count = count + 1;
                end
            end
            
        case 'tca_behav_vis_dist_weighted_temporal'
            tp = datain{1};
            wt = datain{2};
            rgy = datain{3};
            [nc, nw] = size(tp);
            
            count = 1;
%             cols = length(tp{1, 1});
            cols = 4;
            xrg = (-tstep(1) * Fsi: tstep(2) * Fsi) / Fsi;
            skp = 4;
            cnames = {'whnwp', 'nwhnwp', 'whwp', 'nwhwp', 'wh', 'nwh', 'wp', 'nwp'};
            rnames = {'heat w', 'heat nw', 'vf w', 'vf nw'};
            
            figure(gcf)
            clf
            for i = 1: nc
                for j = 1: nw
                    tmp = tp{i, j};
                    for k = 1: cols
                        subplot(nc * nw, cols, (count - 1) * cols + k)
                        [~, id] = max(tmp{k}, [], 2);
                        wtcurr = wt{i, j}(k, :);
                        nt = size(tmp{k}, 2);
                        wtcurr = interp1(linspace(0, nt, length(wtcurr)), wtcurr, linspace(0, nt, nt));
                        wtcurr = wtcurr(id);
%                         wtcurr = wtcurr / sum(wtcurr);
%                         tpp = squeeze(tmp(:, :, k));
%                         [~, idt] = max(tpp, [], 2);
%                         [~, idt] = sort(idt);
%                         imagesc(xrg, 1: size(tmp, 1), tpp(idt, :), [0, 0.2])
                        hold on
                        tt = tmp{k} .* wtcurr';
                        t1 = squeeze(mean(tt, 1));
                        t2 = squeeze(std(tt, [], 1)) / (2 * sqrt(size(tmp{k}, 1)));
                        plot(xrg, squeeze(tt)', 'color', clrs(1, :))
                        plot(xrg, t1', 'color', clrs(1, :), 'linewidth', 2)
                        errorbar(xrg(1: skp: end), t1(:, 1: skp: end), t2(:, 1: skp: end)');
                        xlim([-tstep(1), tstep(2)])
                        ylim(rgy)
                        hold off
                        
                        if k == 1
                            ylabel(rnames{count})
                        end
                        
                        if count == 1
                            title(cnames{k})
                        end
                    end
                    count = count + 1;
                end
            end
            
        case 'tca_avg_temp_vis'
            tp = datain{1};
            se = datain{2};
            if length(datain) < 3
                rgy = [0, 0.2];
            else
                rgy = datain{3};
            end
            [nc, nw] = size(tp);
            
            count = 1;
%             cols = length(tp{1, 1});
            cols = 4;
            xrg1 = (-tstep(1) * Fsi: tstep(2) * Fsi) / Fsi;
%             xrg2 = (-0.5 * Fsi: 1.5 * Fsi) / Fsi;
            xrg2 = (-tstep(1) * Fsi: tstep(2) * Fsi) / Fsi;
            xrg = {xrg1, xrg2};
            skp = 4;
            cnames = {'full whisker', 'no whisker'};
            rnames = {'heat', 'vf'};
            
            figure(gcf)
            clf
            for i = 1: nc
                for j = 1: nw
                    tmp1 = tp{i, j};
                    tmp2 = se{i, j};
                    lgd = [];
                    for k = 1: cols
                        subplot(nc, nw, count)
                        hold on
                        t1 = tmp1(k, round((xrg{i} + tstep(1)) * Fsi + 1));
                        t2 = tmp2(k, round((xrg{i} + tstep(1)) * Fsi + 1));
                        ht = plot(xrg{i}, t1, 'color', clrs(2 * (k - 1) + 1, :), 'linewidth', 2);
                        lgd = [lgd, ht];
                        x = [xrg{i}, fliplr(xrg{i})];
                        y = [t1 + t2, fliplr(t1 - t2)];
                        patch(x, y, clrs(2 * (k - 1) + 1, :), 'edgecolor', 'none')
                        alpha(0.3)
%                         errorbar(xrg(1: skp: end), t1(:, 1: skp: end), t2(:, 1: skp: end)');
                        xlim([xrg{i}(1), xrg{i}(end)])
                        ylim(rgy)
                        hold off
                        
                        if j == 1
                            ylabel(rnames{i})
                        end
                        
                        if i == 1
                            title(cnames{j})
                        end
                    end
                    if i == 1 && j == 1
                        legend(lgd, {'whnwp', 'nwhnwp', 'whwp', 'nwhwp'})
                    end
                    count = count + 1;
                end
            end
            
        case 'tca_behav_subcluster_across_condition'
            tp = datain{1};
            if length(datain) < 2
                rg = [0, 1];
            else
                rg = datain{2};
            end
            [nc, nw] = size(tp);
            
            count = 1;
%             cols = length(tp{1, 1});
            cols = 4;
            nn = 4;
            xrg = (-tstep(1) * Fsi: tstep(2) * Fsi) / Fsi;
            skp = 4;
            lnames = {'whnwp', 'nwhnwp', 'whwp', 'nwhwp', 'wh', 'nwh', 'wp', 'nwp'};
            cnames = {'1st [0-1s]', '2nd [1-2s]', '3rd [2-3s]', '4th [3-5s]'};
            rnames = {'heat w', 'heat nw', 'vf w', 'vf nw'};
            
            figure(gcf)
            clf
            for i = 1: nc
                for j = 1: nw
                    tmp = tp{i, j};
                    for k = 1: nn
                        subplot(nc * nw, nn, (count - 1) * nn + k)
                        p = [];
                        hold on
                        for ii = 1: cols
                            t1 = tmp{ii}(k, :, 1);
                            t2 = tmp{ii}(k, :, 2);
                            x = [xrg, fliplr(xrg)];
                            y = [t1 + t2, fliplr(t1 - t2)];
                            patch(x, y, clrs(mod(ii - 1, 8) + 1, :), 'edgecolor', 'none')
                            alpha(0.3)
                            p(ii) = plot(xrg, t1, 'color', clrs(mod(ii - 1, 8) + 1, :), 'linewidth', 1);
                        end
                        xlim([-tstep(1), tstep(2)])
                        ylim(rg)
                        hold off
                        
                        if k == 1
                            ylabel(rnames{count})
                        end
                        
                        if count == 1
                            title(cnames{k})
                        end
                        
                        if count == 1 && k == 1
                            legend(p, lnames(1: cols))
                        end
                    end
                    count = count + 1;
                end
            end
            
        case 'tca_behav_subcluster_poll_bar_graph'
            tp = datain{1};
            [nc, nw] = size(tp);
            
            count = 1;
%             cols = length(tp{1, 1});
            cols = 4;
            xrg = (-tstep(1) * Fsi: tstep(2) * Fsi) / Fsi;
            skp = 4;
            lnames = {'whnwp', 'nwhnwp', 'whwp', 'nwhwp', 'wh', 'nwh', 'wp', 'nwp'};
            cnames = {'1st [0-1s]', '2nd [1-2s]', '3rd [2-3s]', '4th [3-5s]'};
            rnames = {'heat w', 'heat nw', 'vf w', 'vf nw'};
            
            figure(gcf)
            clf
            for i = 1: nc
                for j = 1: nw
                    tmp = tp{i, j};
                    [nn, nt, ~] = size(tmp{1});
                    subplot(nc, nw, count)
                    hold on
                    
                    t1 = permute(reshape(cell2mat(cellfun(@(x) x(:, :, 1), tmp(1: cols), 'uniformoutput', false)), nn, [], nt), [1, 3, 2]);
                    t1 = squeeze(max(t1, [], 2));
                    barh(-t1(:, 1: 2), 'stacked')
                    barh(t1(:, 3: 4), 'stacked')
                     
                    xlim([-1, 1])
                    ylim([0, nn + 1])
                    xticks([-0.5, 0.5])
                    xticklabels({'no wipe', 'wipe'})
                    if i == 1 && j == 1
                        legend(lnames)
                    end
%                     ylim(rg)
                    hold off
                    title(rnames{count})
                    
%                     if k == 1
%                         ylabel(rnames{count})
%                     end
%                     
%                     if count == 1
%                         title(cnames{k})
%                     end
%                     
%                     if count == 1 && k == 1
%                         legend(p, lnames(1: cols))
%                     end
                    count = count + 1;
                end
            end
            
        case 'tca_model_error'
            dtuse = datain{1};
            a = tca_analysis({dtuse}, 'compute_res_curve');
            a = a{1};
            figure(gcf)
            clf
            hold on
            t1 = [];
            t2 = [];
            for i = 1: 100
                t1(i) = median(a(a(:, i) < 1, i), 1);
                t2(i) = std(a(a(:, i) < 1, i), [], 1) / 2;
            end
            plot(t1)
            errorbar(1: 100, t1, t2, '.', 'markersize', 10)
            xlabel('Component #')
            ylabel('Relative Error')
            title('Model Error Curve')
            
            
            

            
            
            
            
    end
    
end